Magnetic reconnection and stochastic plasmoid chains in high-Lundquist-number 

plasmas 
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A numerical study of magnetic reconnection in the large- Lundquist-number (S), plasmoid- 
dominated regime is carried out for S up to 10 7 . The theoretical model of Uzdensky et al. [Phys. 
Rev. Lett. 105, 235002 (2010)] is confirmed and partially amended. The normalized reconnection 
rate is E e a ~ 0.02 independently of S for S S> 10 4 . The plasmoid flux (vp) and half-width (w x ) 
distribution functions scale as /(^) ~ and f(w x ) ~ w x 2 ■ The joint distribution of $ and w x 
shows that plasmoids populate a triangular region w x > ^ /Bq, where Bq is the reconnecting field. 
It is argued that this feature is due to plasmoid coalescence. Macroscopic "monster" plasmoids with 
w x ~ 10% of the system size are shown to emerge in just a few Alfven times, independently of S, 
suggesting that large disruptive events are an inevitable feature of large-,? reconnection. 

PACS numbers: 52.35.Vd, 94.30.Cp, 96.60.lv, 52.35.Py 



Introduction. It has become clear in recent years that 
resistive magnetic reconnection at asymptotically high 
Lundquist numbers (S) is a temporally and spatially ir- 
regular process, dominated by multiple plasmoids gener- 
ated in unstable current sheets JjJ—lsJ] . The reconnection 
rate in this regime is independent of S provided S > S c 
U, where S c ~ 10 4 H, |I(J is the plasmoid instability 



threshold. Thus, the classic Sweet-Parker (SP) the- 
ory [ll|, [13] is no longer sufficient even for resistive MHD 
reconnection and a new physical paradigm is needed. 

Such a theory was recently attempted by [l3j] (hence- 
forth ULS). The physical picture on which it is based 
is that, as the plasmoid instability [![ proceeds into 
its nonlinear stage, inter-plasmoid current sheets form, 
which arc then subject to the same instability. The re- 
sult is a multiscale plasmoid chain originally envisioned 
by 14|. ULS assume that (i) the current sheets con- 



necting the plasmoids in this chain are typically just 
marginal with respect to the plasmoid instability and so 
their length is ~ L c = (j]/Va)S c , where n is the mag- 
netic diffusivity and Va the Alfven speed based on the 
upstream magnetic field Bq; (ii) the reconnecting field 
is equal to the upstream field Bq for all interplasmoid 
layers and so outflows into all plasmoids are Alfvcnic 
with the same speed Va] and (iii) smaller plasmoids 
do not have time to saturate before they are ejected 
into larger ones (and are promptly merged with them). 
ULS then show that (i) the effective reconnection rate is 

~ —1/2 

E e g = cE c s/B Va ~ S c ~ 0.01; (ii) the plasmoid flux 
(fy) and cross-sheet half- width (w x ) distribution func- 
tions are /(*) ~ E cS BqL^^ 2 and f(w x ) ~ E e sLw~ 2 
(the power laws are the same because it is argued that 
\& ~ w x Bo); and (iii) anomalously large "monster" plas- 
moids occasionally occur, with sizes ~ Sc 1 ^ 4 L ~ 0.1L, 



where L is the system size. Note that diagnosing the 
plasmoid chain in terms of the flux and half-width dis- 
tributions is a natural statistical description for such an 
object 
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lq . Note also that the prediction of mon- 
ster plasmoids is potentially an important one in light 
of the evidence of violent abrupt events associated with 
reconnection sites (e.g., solar flares 17 1 or sawtooth flj|). 

In this Letter, we present a numerical study of resis- 
tive MHD reconnection at the highest currently achiev- 
able Lundquist numbers. Our results confirm the basic 
predictions of the ULS theory, but also reveal that the 
picture is more complex than originally envisioned. 

Numerical setup. We use the same numerical scheme 
as [H| to solve the standard set of compressible visco- 
rcsistivc MHD equations in a 2D box [—L X ,L X ] x 
[—L y ,Ly\. Our setup is designed so that a statistical 
steady state can be reached. Namely, the density, pres- 
sure and the incoming magnetic field are imposed at 
the upstream boundaries (x = ±L X ): p = 1, P = 3 
and By = B Q {l + cos[(7ry/2i + e) 2 ]} /2, where the code 
units are based on Va = Bq/^/Akp = 1 and L = 1. The 
small perturbation e = 0.06L y /L is necessary to break 
the y-symmetry of the numerical set up and thus pre- 
vent the artifical lingering of plasmoids at the center of 
the sheet. Solcnoidality is used to fix B^/x = — B /y at 
the upstream boundary. For the velocity at this bound- 
ary, we set u y /x = 0, whereas u x (x = iL x ) is obtained 
from the frozen-flux condition (the box is wide enough 
that the resistive term is negligible at x = ±L X ). Free- 
outflow boundary conditions are imposed at the down- 
stream (y = ±L y ) boundaries. The initial condition is 
designed to mimic qualitatively a SP-like current sheet. 
This is not, however, a steady-state solution of the resis- 
tive MHD equations, so there is no need to add a pertur- 
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FIG. 1: Top panel: the plasmoid chain for the run with S = 10 6 (16384 2 grid points). Only a fraction of the x-domain is 
shown, —0.03 < x/L < 0.03. The color scale (blue to red) represents B y £ [—1, 1]. Note that the assumption of ULS [Jjl] that 
the reconnecting field is equal to the upstream field Bo all the way to the thinnest of the current sheets appears to hold true. 
Bottom panel: outflow velocity u y (x — 0,y). The outflows into most plasmoids are approximately Alfvenic. 



bation to the initial configuration in order to trigger the 
plasmoid instability. The instability threshold for this 
setup is found to be S c « 1.2 x 10 4 . 

We perform a Lundquist-number scan in the range 
300 < S < 10 7 . In all cases, the viscosity v = 7]. 
Most of our runs are done in a "semi-global" setup with 
L x = 0.3L and L y = 0.5L; the exceptions, for lack 
of sufficient computational resources, are the runs with 
S = 3 x 10 6 (L x = 0.15L^L y = 0.25L) and S = 10 7 
(L x = 0.01L, L y = 0.02L) [251 ] ■ The numerical resolution 
depends on S, ranging up to 16384 2 for S = 10 6 , 3 x 10 6 
and 4096x 8192 for S = 10 7 (for which the box is smaller). 

Reconnection rate. In all our simulations with S > 
S c , the initial SP-like configuration is quickly replaced 
with a plasmoid-dominated current sheet (Fig. [T]) . The 
system then enters a statistical steady state, with mul- 
tiple plasmoids constantly being formed, coalescing and 
being ejected through the outflow boundary. We define 
the effective global reconnection rate in terms the inflow 
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FIG. 2: Reconnection rate E e g [Eq. fTj), squares], and the 
(half) rates of resistive Q v (circles) and viscous Q v (trian- 
gles) heating. Simulations with S > 10 6 last for shorter times 
before they are disrupted by monster ejections and so con- 
verged mean values for heating rates could not be obtained; 
the reconnection rate, calculated at the inflow boundary, did 
not have this problem. 



plasma velocity at the upstream wall: 

where (...) denotes time average. This is plotted in 
Fig. [2] as a function of S. A transition is manifest from 
the SP scaling E cS - S' 1 ' 2 for S < 10 4 - S c to 

— 1/2 

E e g « 0.02 ~ S c , consistent with the ULS predic- 
tion [l3| and previous numerical results 0, 0, • Note 
that this result is now extended to larger values of S than 
ever before. Such an extension is important: as shown 
by ULS, S - Sc /2 - 10 6 is the threshold at which an 
individual plasmoid can saturate faster than it is ejected 
from the global current sheet. This would slow down 
reconnection were it not for plasmoid ejection: smaller 
plasmoids are swallowed by larger ones before they have 
time to saturate. It was assumed by ULS that this coales- 
cence process would operate efficiently — the persistence 
of fast reconnection beyond S ~ 10 6 demonstrated here 
suggests that this assumption is indeed valid. 

Heating rate. The normalized resistive and viscous 
heating rates are Q v = (JJdxdyrij 2 (x,y)/2LyB 2 VA), 
and Qi/ = (ffdxdyvuj 2 (x,y)/2LyVjl) 1 where j z and w z 
are the current and vorticity, respectively. These rates 
are also plotted in Fig. [2] In the fast-reconnection regime, 
Qr) ra Qv ~ 0.008. Since the total Poynting flux into the 
box is (per unit length) P m rj 2E e s ps 0.04 and the ki- 
netic energy influx is small (oc u x E^ s ), the conclusion is 
that ~ 40% of the incoming (magnetic) energy is dissi- 
pated into heat (the rest goes into the reconnected field 
and the kinetic energy of the mass outflows). 

Plasmoid distribution. The plasmoid population is 
naturally characterized by the distribution of fluxes (\E0 
and half- widths (w x ) of individual plasmoids fl3l fl5j ] . 
The distribution functions /("F) and f(w x ) are plotted 
in Fig. [3l The >F _2 and w~ 2 scalings predicted by ULS 
do hold, although the distributions flatten for * / B L 
and w x /L below certain values that decrease at larger 5. 

A more detailed diagnostic is the joint distribution 
function f(^,w x ), which is shown in Fig. |4] and reveals 
a new feature: ULS argued that the plasmoid half-width 
and flux should be related by w x ~ fy/Bo; in fact, 
there is a significant off-diagonal plasmoid population 
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FIG. 3: Plasmoid flux (left) and half-width (right) distribution functions (note that the run with S = 3 x 10 had a shorter 
box by half and so its distributions cut off at smaller flux and width). Dashed lines show the ULS scalings [l3| . 



with w x > ^ I Bq (cf. [Hj]). The presence of these plas- 
moids in the measured distribution can be explained as 
follows. The ULS argument assumed effectively that once 
a smaller plasmoid is ejected into a larger one, it is imme- 
diately and completely absorbed by (i.e., coalesces with) 
the latter and so falls out of the distribution. However, 
in reality, the coalescence between two plasmoids is not 
instantaneous (cf. [l5|]) — and so at any given time, there 
are many plasmoids for which coalescence has started at 
some earlier time and that are in an advanced stage of 
being digested by a bigger plasmoid. Thus, a typical 
plasmoid's life consists of two distinct phases: the ULS 
growth phase (while the plasmoid moves through its host 
current sheet) and the subsequent phase of digestion by 
a bigger plasmoid — this will have an effect on the plas- 
moid distribution. 

We envision the coalescence as a gradual stripping of 
the outer layers of the smaller plasmoid so the magnetic 
field in a semidigestcd plasmoid is B ~ Bqw x /w x q, where 
w x o is the plasmoid's half-width at the beginning of the 
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FIG. 4: Joint distribution of plasmoid flux and half-width for 
S — 10 6 . The solid diagonal line shows the ULS plasmoids 
w x = "if /Bq; the dotted line is the condition |2]). The predator 
and prey plasmoids shown in Fig. [5] are marked by x and +, 
respectively. 



coalescence [26} Its flux is, therefore, ^ ~ Bow x /w x o- 
Since w x < w x q, these plasmoids are off-diagonal: w x > 
^ /Bo. Fig. [5] illustrates the swallowing of a smaller plas- 
moid by a larger one; as shown in Fig. [4j the latter is rel- 
atively close to the diagonal, while the former is strongly 
off-diagonal. 

Let us estimate the widths of the off-diagonal plas- 
moids. These have to be relatively small because larger 
plasmoids take a longer time to be digested and if that 
time exceeds the typical time ta ~ L/Va for both 
predator and prey plasmoids to be ejected from the 
global sheet, then the effect on the measured distri- 
bution is small. The characteristic coalescence time 
is t c \ ~ ^q/cE ~ Bow x o/cE, where is the initial 
flux, cE ~ VaBq max(S' c , Su, 1 ^ 2 ) is the reconnection 
rate, S w ~ VAW x o/r) is the Lundquist number associ- 
ated with the (vertical) current sheet that forms be- 
tween two coalescing plasmoids and we are taking into 
account that reconnection rate is independent of S w for 
S w > S Cl or w x q > L c (length of the longest pos- 
sible plasmoid-stable layer (l3 |: cf. [H}). Therefore, 
t c i/TA ~ Sc /2 {w xQ /L)mm(l, ^w x0 /L c ) < 1, or w xQ < 

L max(5'^ 1 ^ 2 , S~ 1 ^ 3 ), is the condition for scmidigested 
plasmoids to contribute to the off-diagonal part of the 
distribution. Since ^ ~ Bqw x /w x q, this translates into 

w x /L < (*/B £) 1/2 max(5- 1 / 4 , S^ 6 ). (2) 
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FIG. 5: Example of coalescing plasmoids: zoom on the right- 
most part of Fig. [T] Same color scheme is used. Lines of 
constant magnetic flux are also shown. 
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FIG. 6: The half-width of the largest plasmoid in the box vs. 
time. Time is measured from the start of steady-reconnection- 
rate phase. The dotted line is the monster threshold 0.05L. 
Inset: time tm to reach the monster threshold vs. S. 

This indeed appears to capture the maximum of /(\&, w x ) 
rather well (see Fig. 2] for S = 10 6 ; similarly good agree- 
ment was found for other values of S). Since ^ must 
be consistent with w x > ^/Bo, the off-diagonal plas- 
moids only matter if w x /L < max(S , c S* -1 / 3 ). Note 
that the coalescence rate becomes independent of rj for 
S > Sc /2 ~ 10 6 . 

Monster plasmoids. The following argument follows 
ULS in a somewhat expanded and amended form. 
Because the flows carrying both flux and embedded plas- 
moids out of the current sheet are roughly linear, u y ~ 
Vau/L (see Fig. [I}, the ejection time for a plasmoid born 
at some location yo in the sheet is < C j = dy/u y ~ 
TAln(L/yo) (this is true both for the global sheet and 
also for any local one, in which case L would be the typ- 
ical length for the latter). Therefore, plasmoids born 
near the center of the sheet remain in the game logarith- 
mically longer than others. While this only leads to a 
logarithmic correction for their flux, 'F ~ E e gVABot e j ~ 
E c ffBoLhiL/yo, the enhancement of their area is much 
greater. A plasmoid grows in area by absorbing all the 
plasma and smaller plasmoids from roughly up to the 
midpoint of the layer that connects it to its neighbor of 
a similar size. One can then see that the plasmoid area 
A grows according to 



dA 
~dt 



Ay(i) 



d * 
dt.B~^ 



Ay(0)e^ TA E eg VA 



(3) 



where Ay(t) is the (exponentially stretched) half-distance 
to the neighboring plasmoid. Integrating ((3j> up to t = i e j 
gives A(t ej ) ~ E cS Ay{0)L(L/y o - 1). 

If the plasmoid was born away from the center of the 
sheet, yo ~ L, then A ~ E c s Ay(0)L and so w x ~ 
A/wy ~ E c ffL ~ O.Oli. We have estimated the y-extent 
of the plasmoid as w y ~ Ay(0), which does not change 
as long as w x < w y . In contrast, for centrally born plas- 
moids, yo *C L. we have w x ~ E e gL 2 /yo at ejection, 



provided w x < w y ~ Ay(0). If the latter condition is not 
satisfied, i.e., if y < E c rL 2 / Ay(0), the plasmoid will 

w y (which will happen 



be circularized as soon as w x 
before ejection) and so its half- width at ejection will be 

w x ~ 4 1 / 2 L(Ay(0)/y ) 1/2 . Since Ay(0) < y , the max- 

~l/2 

imum half-width achievable is w x . mllx ~ E g L ~ 0.1L. 
This is a nearly macroscopic size — the plasmoids that 
reach it were dubbed monster plasmoids by ULS. Only 
those plasmoids stand a chance of achieving monster sta- 

~ 1/2 

tus that are born at yo < E c ^ L ~ 0.1 L. This must be 
consistent with yo > L c (shorter sheets are stable), which 
implies that monsters will only appear if S > Sc^ 4 ~ 10 5 . 

Fig. [5] shows the half-width of the largest plasmoid, 
w^max in the simulation box vs. time. Exponential 
growth to the monster size is manifest — this is defined 
here, somewhat arbitrarily, as u^max = 0.05L (0.1L is 
never actually reached in our simulations, but it is, of 
course, no more than an order-of-magnitude estimate; 
also our simulation domain is smaller than the system 
size, L y < L). This size is usually achieved by just one 
plasmoid at a time, just before it is ejected, whereupon 
Wz.max dips, then recovers as a new monster emerges, 
and so on. The time £m for a plasmoid system to pro- 
duce and grow a monster can be estimated simply as the 
ejection time for a plasmoid born in the relevant central 

~ —1/2 

part of the sheet: t^j ~ i C j ~ ta In E cS , which amounts 
to a few Alfven times, independent of S — this is borne 
out by the numerical results (Fig. [HI inset). For monster 
plasmoids, ^ < B w x (like for the coalescing ones), so 
they occupy the top right corner of the (^! 7 w x ) plane in 
Fig. [4] (note that the large plasmoid in Fig. [5] is a mon- 
ster in the making) . The probability of finding a monster 
(defined by w x > 0.05L) hovers between 1% and 3%. 

Conclusions. We have found that resistive MHD re- 
connection is fast, its rate cE/BqVa = E e s ~ 0.02, in- 
dependently of S. While a similar conclusion has been 
reported before 0, 0, 041] > our study is the first to probe 
Lundquist numbers significantly exceeding the critical 
threshold of 10 6 [l3| in order to show that plasmoid satu- 
ration does not shut down fast reconnection in the high- 
Lundquist-number, plasmoid-mcdiated regime. It also 
confirms that reconnection occurs via a multiscale plas- 
moid chain Tj^- 15, [l9[, characterized by local Alfvenic 
outflows and many coalescing plasmoids. 

Statistics of this "plasmoid turbulence" are measured 
for the first time in terms of the flux- width joint dis- 
tribution — a natural choice both from the theoretical 
[HI, [IH and observational [l|| perspective. The ULS 
scalings w x ~ ^/B , /(*) ~ *~ 2 , f(w x ) ~ w~ 2 are cor- 
roborated, but we also find a substantial "off-diagonal" 
(w x > *i>/Bo) plasmoid population for w x < E c ffL. The 
excess of plasmoids of relatively large size and small flux 
is explained by considering the coalescence between plas- 
moids. Thus, the full picture of the plasmoid "turbu- 
lence" involves not just multiple reconnection sites along 
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the global layer, but also many transverse layers between 
coalescing plasmoids (these layers can themselves break 
up into plasmoid chains [l^L 

Another large-size low-flux subspecies is the "monster" 
plasmoids, also theoretically anticipated by [li|. They 
are born in the middle tenth of the global layer and grow 
to nearly macroscopic size in just a few Alfven times, in- 
dependently of the Lundquist number. This inevitable 
and relatively frequent nature of what can be very vio- 
lent and disruptive events (ejection of a monster from the 
global layer) is reminiscent of the observed bursty char- 
acter of plasmoid ejections in solar flares [I?], lit! and 
perhaps also of the sawtooth crash in tokamaks [la] . 

These results show that even 2D MHD resistive recon- 
nection contains a wealth of strongly nonlinear, stochas- 
tic behavior — a type of MHD turbulence that is only 
now starting to be studied quantitatively. It is encour- 
aging that the simple phenomenology of the ULS model 
131 ] appears to capture some of the essential properties of 
such systems, but it is also now clear that the full picture 
will require a deeper and more quantitative understand- 
ing of plasmoid coalescence and of the extreme events 
such as the emergence of monsters. Finally, many further 
complications will have to be taken into account before 
idealized models can truly describe the real-world recon- 
nection in its full splendor: t^g., kinetic physics d , l2l"l . |22| . 
background turbulence HH3], 3D effects [Hill). 
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Even in this case, the box is still much wider (in x) than 
the thickness of the SP layer, 8sv/L ~ S~ 1/2 ~3x 10" 4 , 
and much longer (in y) than L c /L ~ S c /S ~ 10 -3 . 
When a typical plasmoid is ejected from its host layer 
into a bigger plasmoid, it is long and thin: w x ~ w v E e f{ 
[l3| . However, immediately upon ejection, it is squashed 
against the bigger plasmoid and becomes circularized 
(w x ~ U) y ), while preserving its flux and area. Fhus, just 
before it starts coalescing with the bigger plasmoid, its 

~ —1/2 

width jumps up by a factor of E cS ~ 10, and its mag- 

~ 1/2 

netic field becomes BoE a g to preserve flux. Fhis immedi- 
ately moves the plasmoid vertically upwards by a factor 

~ — 1/2 

E cS from the ULS diagonal. Our arguments can be 
modified to account for this effect — the result is to raise 
the threshold ([2J upwards by a factor of order unity. 



